####################
####################
####################
####
####
####
####
#### Replication File for 
#### Powell and Grimmer (2016)
#### ``Money in Exile"
####
####
#####################
#####################
#####################


load('ExileData1.RData')


##defining length
len<- length

##creating Figure 1
rep1<- lm(over.post~over.pre + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), data = exile_data)

rep1.coef<- rep1$coef[3] + c(-1.96*sqrt(diag(vcov(rep1)))[3], -1.28*sqrt(diag(vcov(rep1)))[3],0,1.28*sqrt(diag(vcov(rep1)))[3] ,1.96*sqrt(diag(vcov(rep1)))[3])



part1<- lm(I(over.post.party)~I(over.pre.party) + exiled*retire + rel_rank + prev_vote + dist_char + as.factor(congs) + as.factor(comm), data = exile_data)

over.bail<- part1$coef[3] + c(-1.96*sqrt(diag(vcov(part1)))[3],-1.28*sqrt(diag(vcov(part1)))[3], 0,1.28*sqrt(diag(vcov(part1)))[3], 1.96*sqrt(diag(vcov(part1)))[3])




par(las=1)
par(cex.axis=1.5)
par(cex.lab =1.5)
par(mar=c(5, 8, 3, 2))
plot(c(0,1)~c(0,1), pch='', xlab='Effect of Exile, Dollars', ylab='', axes=F, frame.plot=T, ylim=c(-0.25,1.25), xlim=c(-5000, 205000))
axis(2, c(0,1), c('Co-partisan', 'Overall'))
axis(1, c(0,50000,100000, 150000, 200000), c('$0', '$50,000', '$100,000', '$150,000', '$200,000'))
for(z in c(25000, 50000,75000, 100000, 125000, 150000, 175000, 200000)){
	arrows(z, -10, z, 10, lty=2, lwd=1, col=gray(0.8), len=0)
	}
arrows(-1e5, 0, 1e10, 0, lty=2, col=gray(0.8), lwd=1, len=0)	
arrows(-1e5, 1, 1e8, 1, lty=2, col=gray(0.8), lwd=1, len=0)
arrows(0, -1, 0, 10, lty=2,lwd=2, col=gray(0))
arrows(rep1.coef[1], 1, rep1.coef[5], 1, len=0, lwd=2)
arrows(rep1.coef[2], 1, rep1.coef[4], 1, len=0, lwd=4)

points(rep1.coef[3], 1, pch=20, cex=2)



arrows(over.bail[1], 0, over.bail[5], 0, len=0, lwd=2)
arrows(over.bail[2], 0, over.bail[4], 0, len=0, lwd=4)

points(over.bail[3], 0, pch=20, cex=2)


##performing robustness checks

##removing the intelligence and budget members

rep1_drop2<- lm(over.post~over.pre + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), subset=which(comm!= 115 & comm!= 242 ), data = exile_data)

rep1_party_d2<- lm(over.post.party~over.pre.party + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), subset=which(comm!= 115 & comm!= 242 ), data = exile_data)

##removing instances of an inversion

rep1_drop3<- lm(over.post~over.pre + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), subset=which(inv!=1), data = exile_data)

rep1_or_party_d3<- lm(over.post.party~over.pre.party + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), subset=which(inv!=1), data = exile_data)



##adjusting for inflation
rep1_adj<- lm(over.post.adj~over.pre.adj + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), data = exile_data)

rep1_party.adj<- lm(over.post.party.adj~over.pre.party.adj + exiled*retire + rel_rank+ prev_vote + dist_char + as.factor(congs) + as.factor(comm), data = exile_data)

##now replicating Figure 2 

ways.1<- lm(fin.post~fin.pre + exiled*retire + as.factor(congs) + rel_rank + prev_vote + dist_char, subset=which(comm==196), data = exile_data)
nways<- lm(fin.post~fin.pre + exiled*retire + as.factor(congs) + rel_rank + prev_vote + dist_char, subset=which(comm!=196), data = exile_data )




par(las=1)
par(mar=c(5, 7, 3, 2))
par(cex.lab = 1.25)
par(cex.axis = 1.25)
plot(c(0,1), xlab='Effect of Exile, Dollars', ylab='', ylim=c(-0.25,1.25), axes=F, frame.plot=T, pch='', xlim=c(-250000, 50000))
axis(2, c(1,0), c('Ways & \n Means', 'Not Ways\n & Means'))
axis(1, c(-200000, -100000, 0), c('-$200,000', '', '$0'))

for(z in c(-250000, -200000, -150000, -100000, -50000,  50000)){
	arrows(z, -10, z, 10, len=0, lwd=1, lty=2,col=gray(0.8))
	}
arrows(0, -1, 0, 10, lty=2,lwd=2, col=gray(0))	
arrows(-1e7, 1, 1e7,1,  lwd=1, lty=2, col=gray(0.8))
arrows(-1e7, 0, 1e7, 0, lwd=1, lty=2, col=gray(0.8))


wway1<- ways.1$coef[3] + c(-1.96*sqrt(diag(vcov(ways.1)))[3], -1.28*sqrt(diag(vcov(ways.1)))[3], 0, 1.28*sqrt(diag(vcov(ways.1)))[3], 1.96*sqrt(diag(vcov(ways.1)))[3])

nway1<- nways$coef[3] + c(-1.96*sqrt(diag(vcov(nways)))[3], -1.28*sqrt(diag(vcov(nways)))[3], 0, 1.28*sqrt(diag(vcov(nways)))[3], 1.96*sqrt(diag(vcov(nways)))[3])

func<- function(obj, y){
	arrows(obj[1], y, obj[5], y, len=0, lwd=2)
	arrows(obj[2], y, obj[4], y, len=0, lwd=4)
	points(obj[3], y, pch=20, cex = 2)
	}
func(wway1,1)	
func(nway1, 0)
title(main='Finance PACs Decrease to Ways & Means')


energ.1<- lm(ener.post~ener.pre + exiled*retire  + rel_rank + prev_vote + dist_char, subset=which(comm==128), data = exile_data)
nenerg<- lm(ener.post~ener.pre + exiled*retire  + rel_rank + prev_vote + dist_char, subset=which(comm!=128), data = exile_data)

e1<- energ.1$coef[3] + c(-1.96*sqrt(diag(vcov(energ.1)))[3], -1.28*sqrt(diag(vcov(energ.1)))[3], 0, 1.28*sqrt(diag(vcov(energ.1)))[3], 1.96*sqrt(diag(vcov(energ.1)))[3])

ne1<- nenerg$coef[3] + c(-1.96*sqrt(diag(vcov(nenerg)))[3], -1.28*sqrt(diag(vcov(nenerg)))[3], 0, 1.28*sqrt(diag(vcov(nenerg)))[3], 1.96*sqrt(diag(vcov(nenerg)))[3])

par(las=1)
par(mar=c(5, 7.5, 3, 2))
par(cex.lab = 1.25)
par(cex.axis = 1.25)
plot(c(0,1), xlim=c(-200000, 50000), ylim=c(-0.25,1.25), xlab='Effect of Exile, Dollars', ylab='', axes=F, frame.plot=T, pch='')
axis(2, c(1,0), c('Energy & \n Commerce', 'Not Energy\n& Commerce'))
axis(1, c(-200000, -100000, 0), c('-$200,000', '', '$0'))

for(z in c(-250000, -200000, -150000, -100000, -50000,  50000)){
	arrows(z, -10, z, 10, len=0, lwd=1, lty=2,col=gray(0.8))
	}
arrows(0, -1, 0, 10, lty=2,lwd=2, col=gray(0))	
arrows(-1e7, 1, 1e7,1,  lwd=1, lty=2, col=gray(0.8))
arrows(-1e7, 0, 1e7, 0, lwd=1, lty=2, col=gray(0.8))

func(e1, 1)
func(ne1, 0)
title(main='Energy PACs Decrease to Energy & Commerce')



load('SensData.RData')

attach(exile_sens_data)

library(lme4)
##the committeeFactors.RData contains the information to include committee fixed effects.If interested in the origin of any of the data we
##can provide all code to merge.  
load('committeeFactors.RData')

talk_96<- lmer(big.dep~big.lag + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==1996))


talk_96_sub<- lmer(big.dep~big.lag + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==1996 & comm.fact!= 115 & comm.fact!=242))

talk_96_inv<-  lmer(big.dep~big.lag + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==1996 & big.inv==0))


talk_96_adj<-  lmer(big.dep.adj~big.lag.adj + big.exile*big.sens*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==1996 ))





talk_12<- lmer(big.dep~big.lag + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==2012))

talk_sub<- lmer(big.dep~big.lag + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==2012 & comm.fact!= 115 & comm.fact!=242))




talk_inv<-  lmer(big.dep~big.lag + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==2012 & big.inv==0))




talk_adj<-  lmer(big.dep.adj~big.lag.adj + big.exile*big.sens + big.exile*big.retire + big.prev + big.dist + big.rank  + (1| big.indiv) + comm.fact, subset=which(big.year==2012 ))




exile_96<- fixef(talk_96)[3] 
ex_96c<- c(exile_96 - 1.96*sqrt(diag(vcov(talk_96)))[3],  exile_96 + 1.96*sqrt(diag(vcov(talk_96)))[3])

ex_96_2<- fixef(talk_96)[3] + fixef(talk_96)[18]
ses_96<- sqrt(diag(vcov(talk_96))[3] + diag(vcov(talk_96))[18] + 2*vcov(talk_96)[3,18])
ex_96_2c<- c(ex_96_2 -1.96*ses_96, ex_96_2 + 1.96*ses_96)




exile_12<- fixef(talk_12)[3] 
exile_12c<- c(exile_12 - 1.96*sqrt(diag(vcov(talk_12)))[3], exile_12 + 1.96*sqrt(diag(vcov(talk_12)))[3])
ex_12_2<- fixef(talk_12)[3]  + fixef(talk_12)[20] 
ses_12<- sqrt(diag(vcov(talk_12))[3] + diag(vcov(talk_12))[20] + 2*vcov(talk_12)[3,20])

ex_12_2c<- c( ex_12_2 - 1.96*ses_12, ex_12_2 + 1.96*ses_12)




par(mar = c(5, 7, 3, 1))
par(las = 1)
par(cex.axis = 1.5)
par(cex.lab = 1.5)
par(mfrow=c(1,2))
plot(c(0,1)~c(0,1), xlim=c(-2000, 1000), ylim=c(-0.25,1.25), xlab='Effect of Exile', ylab='', frame.plot=T, axes=F, pch='', main = '104th Cong.')
axis(1, c(-2000,  -1000,  0, 1000), c('-$2,000', '-$1,000', '0', '$1,000'))
axis(2, c(0,1), c('Non\nAccess', 'Access'))
arrows( 0, 10, 0, -10,  lty = 2, col='black', len =0 , lwd=2)

for(z in c(-2000, -1750, -1500,-1250, -1000,-750, -500, -250, 250, 500, 750, 1000)){
	arrows(z, -10, z, 10, lty = 2, col=gray(0.8))}
arrows(-1e7, 0, 1e7, 0, len=0, col=gray(0.8), lty=2)
arrows(-1e7, 1, 1e7, 1, len=0, col=gray(0.8), lty=2)

	points(exile_96[1], 0, pch=20, cex=2)
	points(ex_96_2[1], 1, pch=20, cex=2)
	arrows(ex_96c[1], 0, ex_96c[2], 0, len=0, lwd=4)
	arrows(ex_96_2c[1], 1, ex_96_2c[2], 1, len=0, lwd=4)
axis(4, c(0, 1), c('', ''))

par(mar = c(5, 3, 3, 2))

plot(c(0,1)~c(0,1), xlim=c(-3000, 1000), ylim=c(-0.25,1.25), xlab='Effect of Exile', ylab='', frame.plot=T, axes=F, pch='', main = '112th Cong.')
axis(1, c(-3000, -2000,  -1000,  0, 1000), c('-$3,000', '-$2,000', '-$1,000', '0', '$1,000'))
axis(2, c(0,1), c('', '') )#, c('Non\nAccess', 'Access'))
arrows( 0, 10, 0, -10,  lty = 2, col='black', len =0 , lwd=2)

for(z in c(-3000,-2750, -2500, -2250, -2000, -1750, -1500,-1250, -1000,-750, -500, -250, 250, 500, 750, 1000)){
	arrows(z, -10, z, 10, lty = 2, col=gray(0.8))}
arrows(-1e7, 0, 1e7, 0, len=0, col=gray(0.8), lty=2)
arrows(-1e7, 1, 1e7, 1, len=0, col=gray(0.8), lty=2)

	points(exile_12[1], 0, pch=20, cex=2)
	points(ex_12_2[1], 1, pch=20, cex=2)
	arrows(exile_12c[1], 0, exile_12c[2], 0, len=0, lwd=4)
	arrows(ex_12_2c[1], 1, ex_12_2c[2], 1, len=0, lwd=4)



##finally replicating Figure 4.  
load('NextYear.RData')
##again, avoiding pathological problems with a factor in a data frame.  
load('NewIndiv.RData')
dev.off() ##ensuring return to standard settings


new_dep<- next_year[,1]
new_lag<- next_year[,2]
new_dem<- next_year[,3]
new_arrive<- next_year[,4]
new_sens<- next_year[,5]
new_exile<- next_year[,6]
new_indiv<- as.factor(new_indiv)
store_exile<- lmer(I(new_dep - new_lag)~  new_dem + new_arrive + new_dem*new_sens + new_arrive*new_sens + new_exile*new_sens +  (1|new_indiv))

##writing a function that will create the standard errors
stand_err<- function(obj, covs){
	vcvs<- vcov(obj)
	part1<- 0 
	for(z in 1:len(covs)){
		part1<- part1 + diag(vcvs)[covs[z]]
		}
	if(len(covs)>1){	
	for(z in 1:(len(covs)-1)){
		for(y in (z +1):(len(covs))){
			part1<- part1 + 2*vcvs[covs[z], covs[y]]
			}
		}
		}
	return(sqrt(part1))
	
	}



stay_dem<- sum(fixef(store_exile)[c(1, 2)])
stay_dem_se<- stand_err(store_exile, c(1, 2))

stay_dem_sens<- sum(fixef(store_exile)[c(1, 2, 4, 6)])
stay_dem_sens_se<- stand_err(store_exile, c(1, 2, 4, 6))

new_rep<- sum(fixef(store_exile)[c(1,  3)])
new_rep_se<- stand_err(store_exile, c(1, 3))

new_rep_sens<- sum(fixef(store_exile)[c(1, 3,  4,  7)])
new_rep_sens_se<- stand_err(store_exile, c(1, 3, 4, 7))


stay_rep<- sum(fixef(store_exile)[c(1)])
stay_rep_se<- stand_err(store_exile, c(1))

stay_rep_sens<- sum(fixef(store_exile)[c(1, 4)])
stay_rep_sens_se<- stand_err(store_exile, c(1, 4))

exile_rep<- sum(fixef(store_exile)[c(1, 5)])
exile_rep_se<- stand_err(store_exile, c(1, 5))


exile_rep_sens<- sum(fixef(store_exile)[c(1, 4, 5, 8)])
exile_rep_sens_se<- stand_err(store_exile, c(1, 4, 5, 8))


par(las=1)
par(mar=c(5, 8, 3, 2))
par(cex.axis=1.25)
par(cex.lab=1.25)
plot(c(0,1), c(0,1), xlim=c(-2000, 6000), ylim=c(1, 8), axes=F, xlab='Change in PAC Contributions, Industry Level', ylab='', frame.plot=T, pch='')
axis(1, c(-2000, -1000, 0, 1000, 2000, 3000, 4000, 5000, 6000), c('-$2,000', '-$1,000', '$0', '$1,000', '$2,000', '$3,000', '$4,000', '$5,000', '$6,000'))
axis(2, c(1, 2, 3, 4, 5, 6, 7, 8), c('New Minority\nExiled\n Non-Access', 'New Minority\nExiled\nAccess', 'New Minority\nRemaining\nNon-Access', 'New Minority\nRemaining\nAccess', 'New Majority\nRemaining\nNon-Access', 'New Majority\nRemaining\nAccess', 'New Majority\nNew Arrival\nNon-Access', 'New Majority\nNew Arrival\nAccess'))
for(z in c(1, 2, 3, 4, 5, 6, 7, 8)){
	arrows(-1e10, z, 1e10, z, len=0, lty=2, col=gray(0.8))}

for(z in c(-1000, -500, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000)){
	arrows(z,-10, z ,10, len=0, lty=2, col=gray(0.8))
	}
	
arrows(0, -10, 0, 10, len=0, col='black', lwd=2, lty=2)
func.draw<- function(obj1, obj2,y){
	points(obj1, y, pch=20, cex=2)
	arrows(obj1 + 1.96*obj2, y, obj1 - 1.96*obj2, y, len=0, lwd=2)
	arrows(obj1 + 1.28*obj2, y, obj1 - 1.28*obj2, y, len=0, lwd=4)
	}
func.draw(exile_rep, exile_rep_se, 1)	
func.draw(exile_rep_sens, exile_rep_sens_se, 2)

func.draw(stay_dem, stay_dem_se, 3)
func.draw(stay_dem_sens, stay_dem_sens_se, 4)


func.draw(stay_rep, stay_rep_se, 5)
func.draw(stay_rep_sens, stay_rep_sens_se, 6)


func.draw(new_rep, new_rep_se, 7)
func.draw(new_rep_sens, new_rep_sens_se, 8)





#####we also include the CRP data set of regulated industries, to clarify those that we include in our analysis
load('RegulatedIndustry.RData')

